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PROBABILITY  DISTRIBUTION  OF  SPECTRAL  ESTIMATES 
OBTAINED  VIA  OVERLAPPED  FFT  PROCESSING 
OF  WINDOWED  DATA 

INTRODUCTION 

Estimation  of  the  autospectrum  of  a stationary  random  process  by 
means  of  overlapped  FFT  processing  of  windowed  data  (the  so-called 
direct  method)  is  a popular  and  efficient  method,  especially  for  data 
with  pure  tones  present.  Stable  spectral  estimates,  as  measured  by 
the  equivalent  degrees  of  freedom  of  the  spectral  estimate,  result 
when  the  product  of  the  available  record  length  and  the  desired 
frequency  resolution  (the  time -bandwidth  product)  is  large  in  com- 
parison with  unity.  (See,  for  example,  references  1 and  2 and  the 
references  listed  therein.) 

The  equivalent  degrees  of  freedom  of  the  spectral  estimate  is  an 
incomplete  probabilistic  descriptor,  because  it  depends  only  on  the 
mean  and  variance  of  the  random  variable,  In  this  report,  we  address 
the  problem  of  obtaining  the  characteristic  function  of  the  spectral 
estimate  with  overlap  processing, of  a signal  tone  present  in  Gaussian 
noise,  and  thence  the  cumulative  probability  distribution  (perhaps  by 
numerical  means  as  given  in  references  3 and  4).  For  the  case  of 
signal-absent  also,  we  will  compare  the  exact  probability  distribution 
with  an  approximate  distribution  that  uses  only  the  first  two  moments 
of  the  spectral  estimate,  to  see  when  the  approximate  distribution  can 
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be  used  as  a valid  probabilistic  description.  Some  related  work  is 
available  in  reference  5 and  the  papers  cited  therein, 

A discussion  of  the  relative  stability  of  the  spectral  estimates 
with  signal  tones  present,  and  of  a cross-spectral  estimate, completes 
the  presentation. 

CHARACTERISTIC  FUNCTION  FOR  SIGNAL  PLUS  NOISE 

The  method  and  conditions  of  processing  arc  described  fully  in 
reference  1 and,  for  sake  of  brevity,  will  not  be  repeated  here.  The 
,power  spectral  estimate  at  analysis  frequency,  f,  is  given  by 
(reference  1,  pp.  2-4) 

P 

5(f)  4 V |yt)|!,  (1) 

p=1 

where  P is  the  total  number  of  weighted  data  segments.  Here* 


Y (f)  = 
P 


/ 


dt  exp(-i2irft)  x(t) 


w t - I L - (p-l)S  , 


(2) 


where  x(t)  is  the  available  (complex)  data  process,  w(t)  is  the  data 
window  of  length  L,  and  S is  the  amount  of  shift  each  adjacent  data 
window  undergoes.  The  fractional  overlap  is  therefore  1 - S/L. 


* Integra  Is  without  limits  are  over  the  range  of  non-zero  integrand 

Best  Available  Cor' 
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If  we  let  x(t)  be  composed  of  a pure  signal  tone* 

s(t)  = A exp(i2Ttfot+i0)  (3) 


and  zero-mean  Gaussian  noise  n(t) , (2)  can  be  expressed  as 

Y = Y + Y , 

P PS  pn 

where  the  variable  f is  suppressed  for  notational  convenience  and 
complex  (non-random)  constant 

Yps  = A W(f-fQ)  exp[i0-i2ir(f-fo)(I  L + (p-l)sj  , 


(4) 


(5) 


where 


f 


W(f)  = /dt  exp(-i2irft)  w (t ) i 


(ft) 


| W(f ) | 2 is  called  the  spectral  window  (see  equation  (5) , reference  1), 
and  has  analysis  bandwidth  B.  Now  if  analysis  frequency,  f,  is  not. 
within  a bandwidth,  B,  of  tone  frequency,  f0,  (5)  will  be  zero;  there- 
fore,wc  limit  consideration  to  |f-fQ|<B.  The  remaining  term  in  (4), 


pn 


= exp(-i2irft)  n(t)  w£t  -if-  (p-l)sj 


(7) 


is  complex  Gaussian  since  n(t)  is  Gaussian, 

Substituting  (4)  in  (1),  the  spectral  estimate  is  given  by 

I> 


(8) 


*The  generalization  to  several  separated  tones  will  be  obvious. 
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where  {Y  } are  complex  constants  arid  {Y_„}  are  complex  correlated 
ps  pn 

Gaussian  zero-mean  random  variables, and  the  correlation  is  iuc  to  the 
overlapped  processing. 

In  appendix  A,  the  characteristic  function  of  forms  like  (8)  is 
evaluated;  it  specializes  here  to  the  form 


c(0  = IX 

p*i 


( 1 -iX  f/P)  exp 


, , i ly  |2A  z/p 

I 1 p1  p- 


1-iXpC/P 


(9) 


where  {A^}  are  the  eigenvalues  of  I’  x I1  matrix 


and 


k = C:  Jy  y *| 

) Pn  qn| 

U J 1 I 


y Q K m 


(10) 


HD 


where  Q is  the  normalized  modal  matrix  of  K,  and 

m - TY  •••*  1T 

L Is  PsJ  . (lJ) 

The  evaluation  of  K in  (10)  is  considered  in  appendix  B It  is 

given  by 

K » Ik  ] - g (f)4>  (o)  r , (15) 

L q-pj  n w 

where  Gfi(f)  is  the  noise  spectral  density  at  analysis  frequency,  f; 


■/ 


(T)  = / dt  w(t)w*('t-l  J ; 


( N) 


'1 
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and 


[Vp] 


1 r r . . .r  , 

1 2 P-1 


where 


r ■ 


1-P 


<j>  (mS) 
w 

V°> 


(15) 


(16) 


A Fourier  transformation  of  (9)  would  yield  the  probability 
density  function  of  the  spectral  estimate  (8),  for  a tone  present. 
This  would  have  to  be  done  numerically,  but  has  not  been  pursued  here. 


MEAN  AND  VARIANCE  FOR  SIGNAL  PLUS  NOISE 
By  means  of  (A-16),  the  mean  and  variance*  of  spectral  estimate, 

A 

G(f),  in  (8)  can  be  expressed  as 

Mean  {8(f)}  = Kq  ♦ ImJ2  , <17) 

k=l 

(>-T!)|Kkl1  * i m H K n,  , <18’ 

• k=l -P 

in  terms  of  the  quantities  in  (12)  and  (13).  Employing  the  explicit 
relationships  in  (12)  and  (13),  there  follows 

Mean  {8(f)}  = Gn(f)4>w(0)  ♦ A2|W(f-fo)|2  , (19) 

*More  generally,  the  cumulants  are  given  by  (A-7) . 
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V.,r  {6(f)}  - lrk|! 

k=  1 - P 
P-1 

♦ 2A*  | W(f-fQ)  | aGn(f)  »w(0)p-  ^ /l~  NV  cxp(ik2ir(f-f0)s)  , (20) 

k=l-P  ' ' ' 


where  we  have  employed  (15)  and  (5) , 

At  this  point,  it  is  convenient  to  define  the  output  signal  power 

of  a window  filter  with  transfer  function,  W,  centered  at  f as 

Q (f)  » A2 | W(f- f ) |2  , (21) 

S 0 

and  the  output  noise  power  of  the  same  filter  as 


Vf>  - 


°fiv 


Gn(p)|W(p-f)|2  a 


dp|w(u-f)|2  = 


(;„(f)t(°) 


(22) 


Then  (.19)  and  (20)  take  the  forms 

Mean  |(1(f)J  = Qn(f)  + Qs(0 

k=l-P 


(25) 


+ 2Qs(f)Qn(ijl  T”)  ^ exptlk2l,(f_f°)S)‘  (24) 

k=l  -P 

Prom  (24),  we  see  that  the  presence  of  signal  (A  t 0)  always  increases 
the  absolute  level  of  the  variance  of  the  spectral  estimate  over  that 
for  noise-alone.  If  the  noise  is  absent,  the  variance  of  the  estimate 
is  zero.  If  the  signal  is  absent,  the  equivalent  degrees  of  freedom, 
defined  as 


(> 
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EDF 

n 


2 (Mean) 2 
Variance 


2 


i £(t)v 

k“l-P 


(25) 


is  identical  ^to  equation  (8),  reference  1,  as  it  should  be. 
On  the  other  hand,  for  Qs(f)  » Qn(f), 


_ 2 (Mean) 2 ^ 

EDF  = - — 3 

s Variance 


Qs(f) 


Qnff)  F 2Z  (I_  ^)rk  exP<ik2"tf-fo)S) 


. (26) 


k=l-P 


When  a strong  signal  is  present,  EDFS  is  larger  than  EDFn  by  approxi- 
mately the  ratio  — Qs(f)/Qn(f)  >>  1,  since  the  ratio  of  sums  in  (25) 
and  (26)  is  approximately  unity  for  f a fQ  and  reasonable  overlaps 
(see  (27)  below) , That  is,  the  relative  fluctuation  in  the  spectral 
estimate  is  reduced  by  the  addition  of  signal,  even  though  the  absolute 
variance  increases, 

For  Hanning  weighting  and  50%  overlap  (S  = I./2),  we  find  rQ  ■ 1, 
r = 1/6,  r.  * 0 for  kit  2.  Then  the  two  sums  in  (25)  and  (26)  take  on 

+ 1 K 

the  values 


1 + u-phT  ’ 1 + OpOy cos  [ 2n(f-f0)s]  , 


(27) 


respectively,  The  former  value  is  slightly  larger  than  unity,  whbreas 
the  latter  value  varies  between  approximately  2/3  and  4/3,  depending  on 
the  exact  location  of  the  signal  tone  frequency,  fQ,  with  respect  to 
the  analysis  frequency,  f.  For  an  FFT  approach,  at  least  one  bin  has 
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its  frequency  location,  f,  such  that  |f-fQ|5  (21,)"*;  thus,  at  least 
one  frequency  bin  is  located  such  that  the  latter  value  in  (27)  is 
larger  than  unity. 

Figure  1A  represents  the  power  spectral  estimate,  (1),  plotted  on 
a linear  scale  proportional  to  watts.  The  "ribbon  width"  in  the  region 
of  noise-alone  is  denoted  by  a.  The  amount  of  fluctuation  of  the 
estimate  at  f is  denoted  by  b and  is  larger  than  a.  (The  quantity  b 
is  obsorvable  only  by  rerunning  the  spectral  estimation  procedure  for 
.independent  noise  segments.) 


* A 

G(f)  10  log  G(t) 


Figure  1. 


Spectral  Iistimates  for  Signal  Plus  Noise 


If,  instead,  the  power  spectral  estimate  is  plotted  on  a dll  scale 
(see  figure  IB),  the  noise-alone  ribbon  width,  c,  is  larger  than  the 
fluctuation,  d,  of  the  estimate  at  f . The  mathematics  behind  this 

i 

8 j 

j 

\ 
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conclusion  follows.  Let  the  spectral  estimate  at  frequency  f be 
expressed  as 

5(f)  = ni'-x,  C28) 

where  m is  non-random  and  x has  zero-mean  and  variance  a2,  Then 

<Tb  = 10  log  G (f)  = 10  log  m + 10  log(l+i).  (29) 

Now  suppose  that  a/m«l,  which  could  be  realized  by  means  of  a large 

number  of  pieces,  P,  or  a high  signal  to  noise  ratio;  then 

cfe  a 10  log  m + £ . (30) 

in  10  m 

The  last  term  in  (30)  is  proportional  to  the  relative  stability  of  the 
spectral  estimate  (28);  in  fact 

Var  a (liTo)  ’ (31) 

which  can  be  made  arbitrarily  small.  Thus  a plot  like  figure  1 is 
easily  achievable  and  should  be  anticipated  for  a pure  tone  in  Gaussian 
noise. 

PROBABILITY  DISTRIBUTION  FOR  NOISE-ALONE 

A 

For  noise-alone,  the  mean  and  variance  of  spectral  estimate,  G(f), 
are  available  from  (19),  (20),  and  (16)  as 


Mean  {5(f)}  = Gn(f)^(0), 
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which  agree  with  equations  (5)  and  (6),  reference  1,  respectively. 
More  generally,  the  characteristic  function  follows  from  (9)  as 

P 

p»i 


(S3) 


Now  lot  us  define  a normalized  random  variable 

0(f) 


g 


(54) 


(1  (f)4>  (0)  ’ 

n 'Twv 


notice  that  the  scale  factor  is  independent  of  P and  the  amount  of 
overlap.  Thus  the  mean  li{&)  ~ 1,  and  tin;  characteristic  function  of 

g is 


v°  - 


P-1 


-1 


(35) 


where  are  the  eigenvalues  of  matrix  R in  (15).  Then  by  a 

partial  fraction  expansion,  the  probability  that  random  variable  g 
remains  below  a threshold  value, v,  is  found  to  he 

1’rob  (g<v)  = 1-N  B exp  I-  -'~T~  V v>()' 

^ V >>  ) 

r 1 1 M 

B 


k=  l 


(3ft) 


whe  re 


[‘•a1- 


isksf, 


(37) 


rr  jc  - 


p-i 

pa 
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We  have  assumed  all  the  eigenvalues  of  R to  be  unequal;  this  is  the 
case  if  the  overlap  is  greater  than  0,  which  is  the  case  of  most 
practical  interest.  The  eigenvalues  are  all  non-i  ative  since  R is 
a non-negative  definite  matrix  (see  appendix  B) . 

Equation  (36)  is  an  exact  expression  for  the  cumulative  probability 
distribution  in  terms  of  the  eigenvalues  of  matrix  R.  If  we  consider 
another  random  variable,  t,  with  the  same  mean  and  variance  as  g,  a 
candidate  approximate  characteristic  function  is  (guided  by  form  (35)) 

Ct(5)»  (l-i?/b)'b  , (38) 


where,  in  order  to  maintain  the  same  variance,  we  choose 

P P P-1 

>2 


"..ME  (-^) 

p=l  p,q=l  k=l-P 


(39) 


Equation  (8),  reference  1,  shows  that  b = K/2,  i.e.,  half  of  the 
equivalent  degrees  of  freedom,  Then  the  approximate  probability 
density  function  is 


P(t) 


1 


r(b) 


, b b-1  -bt 
b t e , t>0, 


(40) 


and  the  approximate  cumulative  probability  distribution  is  (equations 
6.S.2  and  6.5.12,  reference  6): 


/ 

0 


dt  p(t)  = 


T(b+1) 


(bv)be“bv  jP  (l;l+b;bv),  v>0.  (41) 
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(A  further  simpler  approximation,  not  pursued  here,  would  he  to  sot 
= integer  part  of  b,  b-,  = bj  + 1,  and  bracket  the  results  above  by 
two  simpler  sums.) 

We  shall  now  make  quantitative  comparisons  between  exact  result 
(36)  and  approximation  (41)  which  has  the  same  mean  and  variance.  The 
question  is:  is  b in  (39)  and  (41)  a sufficient  statistic  to  accurately 
quantitatively  describe  the  exact  cumulative  probability  distribution 
(36),  for  representative  data  windows,  overlap,  number  of  pieces,  and 
time-bandwidth  products,  over  the  range  of  probabilities  of  interest 
to  most  users?  rf  so,  then  attention  can  be  confined  to  the  equivalent 
degrees  of  freedom  and  its  maximization  alone,  as  was  done  in  reference 
1;  this  simplification  would  be  most  worthwhile  and  of  obvious  impor- 
tance . 

The  actual  numerical  computation  of  the  cumulative  probability 
distribution  Prob(2'  v) , is  considered  in  appendix  C . In  figure  2,  the 
exact  cumulative  probability  distribution  for  Hanning  windowing  is 
presented  for  time-bandwidth  product'  BT  = B,  16,  32,  64  , where  T is 
the  available  record  length  and  15  is  the  desired  resolution  bandwidth. 
In  each  plot,  the  overlap  is  varied  from  0 to  approximately  7596,  and 
the  distribution  plotted  on  a normal  probability  ordinate  covering  the 
range  (.0001,  .9999).  The  fact  that  the  curves  are  not  straight  lines 
over  this  range  means  that  a (iaussian  approximation  to  the  power 
spectral  estimate  would  not  suffice.  However,  the  (iaussian  approxima- 
tion would  be  a fairly  good  one  for  larger  BT  and  1’  (see  figure  21),  for 
example)  . 
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The  fact  that  the  curves  in  figure  2 are* virtually  identical  for 
overlaps  greater  than  50%  means  that  there  is  little  point  in  choosing 
overlaps  greater  than  this  amount.  This  confirms  the  choices  of  over- 
lap made  in  reference  1,  where  attention  was  confined  to  the  equivalent 
degrees  of  freedom.  The  ideal  distribution  would  be  a vertical  line 
at  v = 1;  the  closeness  of  these  curves  to  the  ideal  is  a measure  of 
the  spread  of  the  spectral  estimate. 

The  corresponding  results  for  the  approximation  (41)  are  presented 
in  figure  3,  The  curves  are  virtually  identical  to  those  of  figure  2 
over  the  complete  range  of  probabilities  considered,  for  various  values 
of  BT  and  overlap. 

For  a cubic  window,  the  exact  results  and  the  approximation  are 
given  in  figures  4 and  5,  respectively,  The  conclusions  are  identical 
to  those  made  for  the  Hanning  window. 
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FLUCTUATIONS  OF  CROSS  SPHCTRAI,  FSTI MATH 


This  topic  Is  not  directly  related  to  the  earlier  ,'iatoriul  on 
autospectral  estimation;  however,  it  is  an  important  observation  and 
merits  a comment.  For  two  uncorrelated  processes,  x and  y,  the  cross 

A 

spectrum  (!  (f)  ~ 0.  However,  the  cross  spectral  cst  iinate , C ( f)  , 

xy  xy 

satisfies  the  equations  (reference  2): 

k {igifi} » o. 


f:  { | c (fj|2}=|n  (fje;  (f)  s 2a2, 

( 1 xy  1 J K XX  yy 


where  K is  the  equivalent  decrees  of  freedom.  Now,  if  K»l,  1!  (fj  is 

xy 

approximately  complex  Gaussian.  Therefore,  if  we  define  the  amplitude 


estimate 


A S |(i  if) 
' xy 


it:  has  probability  density  function 


Then  the  mean  of  A Is 


(-  -4 

\ la' 


W #<-(: 


it  (lxx  1 * ^'yy  ^ * * 


which  is  it  rather  slow  decay  with  K.  Then  the  ratio  of  the  mean 
amplitude,  (45),  to  the  square  root  of  the  product  of  the  auto-spectra 
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[GxxCfJGyyff)]^ 


1.253  . 

kH 


(46) 


If,  for  example,  X = 32,  this  ratio  is  ,222  which  is  -6,55  dB;  this  is 
not  very  far  down  relative  to  unity  coherence,  though  the  two  processes 
are  uncorrelated. 

Also, 


Var  {A}  = 


g^-^Gxx(f)Gyy(--  , 


(47) 


and,  thereforo, 

Standard  deviation  JaI  „ |±2V*  = 0.52,  f4, 

Mean" '{A}  \ H 1 

independent  of  X (or  P) . So  for  a zero  cross -spectrum  value,  A <» 

|G  (f ) | will  always  have  the  same  amount  of  relative  variation, 
xy 

regardless  of  the  number  of  pieces  P (for  large  P) ; thus,  on  a dB 
scale,  the  "ribbon  width"  of  the  cross-spectral  estimate  is  indepen- 
dent of  P,  when  the  two  processes  are  uncorrelated. 
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DISCUSSION 

An  exact  expression  for  the  characteristic  function  of  the  power 
spectral  estimate  of  a pure  tone  in  Gaussian  noise  has  been  attained, 
and  then  specialized  to  noise-alone,  In  the  noise-alone  case,  a 
numerical  computation  of  the  cumulative  distribution  function  has  been 
conducted.  Comparison  of  the  latter  with  an  approximation  utilizing 
only  the  mean  and  variance  shows  excellent  agreement  over  a wide  range 
of  probabilities,  regardless  of  the  exact  window,  overlap,  or  the  time- 
bandwidth  product.  This  means  that  concentration  on  the  equivalent 
degrees  of  freedom,  particularly  on  its  maximization,  is  sufficient 
for  a probabilistic  description  of  the  auto-spectral  estimate. 

Maximizing  the  equivalent  degrees  of  freedom  results  in  a narrower 
probability  density  function,  as  witnessed  by  the  increased  steepness 
of  the  cumulative  probability  distributions  presented, 

An  entirely  different  method  of  auto-  and  cross -spectral  estimation 
has  been  presented  in  references  7 and  8,  and  is  mentioned  here  as  a 
viable,  attractive  alternative,  particularly  for  short  data  segments. 
Since  only  a few  parameters  are  estimated,  the  estimates  are  potentially 
more  stable,  whereas  the  technique  considered  here  (and  in  reference  1) 
assigns  independent  degrees  of  freedom  to  each  and  every  frequency 
cell  of  interest  and,  therefore,  requires  the  estimation  of  many  more 
parameters. 
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Appendix  A 

DERIVATION  OF  CHARACTERISTIC  FUNCTION 

The  first  half  of  appendix  C of  reference  9 considers  the  Hermitian 

form  ,, 

F =.XH®X,  (A- 1) 

with  mean  and  covariance  of  the  complex  random  variable  matrix  X, 

e)X|  = m , Covjx}  = e j(X-m)(X-m)H|  a K , (A-2) 

where  matrix  X is  P x 1,  and  matrix  B is  Hermitian  and  P x P.  Defining 
P x P matrix 

A " K ^BK^  , CA-S) 

with  corresponding  normalized  modal  matrix  Q and  (diagonal)  eigenvalue 
matrix  we  can  express  (A-l)  as 

P 

Ak|vJ<  , 

k=l 

where  matrix  V Is  P x 1 with  mean  and  covariance 

Mvf  =QHK'5im  5 n,  Cov  |v[  =1  . 

Then  a slight  generalization*  of  the  second  half  of  appendix  C 
of  reference  9 (sec  also  reference  10)  yields  the  characteristic 
function  of  random  variable  F in  (A-4)  as 

*We  must  also  have  E |(X-m)(X-m)j>  = O,  in  addition  to  (A-2). 


(A-4) 


CA-S) 


A-l 
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cm  = 


n 

k=l 


'>  |uj\£ 

d-iXk5)_'  oxpl K— - - 


(A-(>) 


whore  {A^}  and  {]j^j  are  the  elements  of  matrices  \ and  pt  . The  cumulants 
of  F follow  easily  from  (A-6)  as 


c = (n-1) 
n 


l> 


( A - 7 ) 


k=  1 


In  particular,  the  first  two  cumulants  arc 

P 


Xk 


Mean  |f|  = c^ 

,,  k = l 

V»r|F|  . c2  tl*2|pkl*) 


(A-Hj 


For  the  case  of  zero-mean  variables,  L.c.,  itl  = O,  (A-5)  yields 
p = O , and  the  characteristic  function  becomes 

for  zero-mean  variables.  (A-9) 


C(Q  = 

k=l 

The  cumulants  are  then 


sf 


Cl-iAkO 


Cn  . 


for  zero-moan  variables. 


(A- 10) 


k=l 


i, 

'2 


. 1 

pt  is  not  necessary  to  evaluate  K for  eigenvalue  purposes  alone, 
because  the  eigenvalues  { A ^ } of  matrix  Adcfincd  in  (A-3)  arc  the  saun- 
as the  eigenvalues  of  KBor  BK .) 
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As  a specific  application  of  the  general  results  above,  we  consider 


B = I 


[”l  ' 


r K 


[vj 


(A-ll) 


Then  from  (A-3) , we  see  that  A = K . In  order  to  evaluate  (A-8) , we 
notice  that 

£ ^ - £ y ■ • 

k=l  p=l 


k=l 


m 


HK~ij  A K m=mH  m 


(A- 13) 


I’  P P-1 


yv  a a = y 


|K  |2 
p-q 


j>i*  • 

V CP-UI)|KkP  . ca-hj 


k=i  p>qBl 


p.q=i 


k=l-P 


r 


k=l 


V2|u.  I2  = n HXXH  = m 11  k '5  QXXQK 

k k 


m K AQ\Q  hK  ~h  m = mH K 1 KKK* m =mHKw,  (A-iS) 


Then  (A-8)  yields 


Mean  )p|  = = PK  |nijz  ’ 

P-1  k'“'1 

Var  {F}  = c2  (P-jk|)|Kk|z  + 2 ltlH 


Km 


(A— 16) 


k = l-P 


The  specialization  to  zero-mean  variables  is  obtained  by  dropping  the 
last  terms  in  (A-16), 
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Appendix  B 


DERIVATION  OF  COVARIANCE  MATRIX 


We  are  interested  in  deriving  the  two  averages 

EiY  Y and  hiY  Y l (B-l) 

\ pn  qnj  ) pn  qn/ 

because  they  are  needed  for  appendix  A and  to  see  if  the  conditions 
required  there  are  satisfied,  We  have,  from  (7), 


Ii|Y  Y *j,  -j'J dt  du  exp(-i2fff(t-u))  E|n(t)n*(u)|  w|\-“L- (p-l)sj 
w*  ju-|l,-(q-l)sj  . 


(B-2) 


Letting  the  noise  correlation  in  (B-2)  bo  denoted  by  Rn(t-u),  and  its 
spectrum  by  Gn>  (B-2)  becomes 


E jYpnYq*}  = J 'dp  Gn(M) cxp(i2ir(n-f)t)  w |t-~-L- (p-J)sj . 

j/du  exp(i2ir (p-f)u)  w ju-iL- (q-l)sj  j 


■/ 


du  Cn(y) I W(f-M)  I 2 


expj^i2ir (f-p)  (q-p)sj 


(B-3) 


This  quantity  is  a function  only  of  the  difference  of  indices  q and  p. 

If  spectral  window  |W|2  is  narrower  than  the  detail  in  noise 
spoetrum  G[(  in  the  neighborhood  of  analysis  frequency  f,  (B-3)  simplifies 
to 


B-l 


E {Vqn}  3 Gn(f)/dM|l»(M)|*  exp[i2ir(f-M)(q-p)s] 


= Gn(f)<|)w((q-p)S), 

where 

♦W(T)  = /"dt  W(t)W*(t-T) 

is  the  autocorrelation  function  of  data  window  w. 
Now  let 

VraS) 

I^ToT  “ r,n 

Then 

oiy  Y *1  = (1  (f)<J>  (0)  r 

( pn  qnj  n ' rw  q-p 

and  from  (10), 


where  P x P matrix 


K=  Gn(f)<(.w(0)Pt  , 


(B-4) 


(B-5) 


(B-G) 


(B-7) 


(B-8) 


(B-fl) 


is  llerinitian,  Toeplitz,  and  non-negative  definite.*  For  real  weighting 
w,  R is  a real  symmetric  Toeplitz  matrix,  'i’he  matrix  in  (B-8)  is  the 
one  required  in  (A-ll)  and  (10). 

*Thi.s  property  is  easily  proven  by  use  of  definitions  (B-5)  and 
(B-6) . 
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The  second  quantity  we  desire  is,  by  use  of  (7), 

K |YpnYqn|  dt  du  exp(-i2Trf(t+u))  P,|n(t)n(u)|  w^t-^JL- (p-l)sj  . 

w JV^L-Cq-nsJ  . (B'1()) 

Letting  the  noise  correlation  i.n  (B-10)  be  denoted  by(Rn(t-u),  and  its 
spectrum  (Fourier  transform)  by  (J n » (B-10)  becomes 

l'{YpnYqn}  =/d^n(./dt  exp(i2ir(u-f)t)  w[t-±t-(p-l)s]  . 

y*du  cxp(-i27T(p+f)t)  W (u-~L- (q-l)sj 
sj  dy  (J^  (p)W(f-U)W(f+p)  exp£-i2irf  (L-2S+pS+qS)  - j 2'irp (q-p)sj  _ (B-ll) 


If  analysis  frequency  f is  greater  than  the  bandwidth  B of  spectral 
window  W,  then  W(f-p)  and  W(f+q)  do  not  overlap  on  the  p-axis,  and 
(B-ll)  yields 

D/Y  Y \ a 0 if  f>B.  (B-12) 

\ Pn  qnf 

Thus,  the  property  desired  in  appendix  A (footnote  to  equation  (A-6)) 
holds  true  if  f>B. 
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Appendix  C 


NUMERICAL  COMPUTATION  OF  CUMULATIVE  PROBABILITY  DISTRIBUTION 


The  numerical  computation  of  the  cumulative  probability  distribution 
Prob(g<v)  is  not  accomplished  here  directly  via  the  sum  in  (36).  The 
reason  is  that,  for  largo  P,  (36)  is  an  alternating  stun  of  terms  of 
large  magnitude,  and  accuracy  is  lost  in  the  final  result.  Instead, 
the  methods  in  references  3 and  4 are  utilized  on  characteristic 
function  (35):  for  a random  variable  limited  to  positive  values,  the 

cumulative  probability  distribution  can  be  expressed  as  (reference  4) 


OO 


*!«) 

5 


cxp(-i£v) 


, v>0, 


(C-l) 


where  f. (£)  is  the  imaginary  purt  of  the  characteristic  function  f(£). 
We  have  (£)/£*•* E Igl  = 1 as  ( + 0.  We  approximate  (C-l)  according  to 

P(v)  s 1 - — Rej  J 7 d£  exp(-.i£v)  j , (^-2) 


and  then  sample  and  approximate  this  expression  according  to 

2 ( fi  ^ r \ 

P (n  Av)  a J - ^ RejAg^  ^ wr  cxpl-ikAf;  n Avj  > , 

k-o  k 

where  L A£  = £2,  and  <w^>  arc  Trapezoidal  weights  of  integration. 


(C-3) 

Wo 


choose  sampling  increment 

A£  = 

where  N is  chosen  large  enough  that 


2tt 

Fav  ' 

t'i(?)/£  does  not  change  much 


(C  - 4 ) 
i n 


C-l 
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Af).  Then 


' L 


P(n  Av)  - 1 - 


\ Re  jA^  Wk  fi-kAf~  e*p[‘i2lTkn/N]} 


k=0 


N-l 


= 1 


~ AC  Re  | gk  exp  i 2iTkn/N J | 


(C-5) 


where 


k=0 


J 


8 


z 


fi((k+jN)  A?) 

k jLmi  k+jN  (k+jN)  AT “ 
j=0 


C£k<N-l,  ~~-l 


(C-h) 


liquation  (C-5)  is  an  N-point  PIT;  therefore,  we  choose  N as  a power  of  2 
for  speed  purposes. 

The  only  remaining  question  is  the  choice  of  limit  £2  in  (C-2). 
Prom  (35) , we  know  that 

|f.(C)|  < |f(Q|  i 1 (C-7) 

P 

n {*•  rP5} 

p=i 

where  r_  = Therefore 

P p 


iMOl  < 


1 


(<:-«) 
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where  IP  can  be  1 or  2 or  ...  or  P.  Therefore  the  error,  E,  in  using 
(C-2)  rather  than  (C-l)  is  bounded  according  to 


This  equation  can  be  solved  for 


2 

7T 


nw 


-p+ip-i 

S-2 

P+1-IP 


(C-9) 


5*  - •* 


TTE 


1 

2 P+l-IP 

P 

.nwi^  . (c-io) 

p=ip 

with  the  guarantee  that  the  error  will  be  less  than  E if  we  choose  £2 

according  to  (C- 10) . Since  IP  is  not  unique,  we  choose  £?.  to  be  the 

minimum  value  over  the  range  of  IP=1,  2,...,  P,  for  then  the  integration 

range  in  (C-2)  can  be  Kept  to  a minimum. 

In  summary,  we: 

spec) fy  Av , E , P,  BT 

compute  and 

choose  N = 1024  (say) 


2*ir 

compute 

compute  L = £?./A£; 
let  .1  =>  (b+D/N 
compute  (C-6) 

compute  ITT  <jgkj  and  printout  (C-5) 

choose  N = 2048,  go  back  to  step  4,  and  observe  change  in  printout. 
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A program  for  this  procedure  fox  the  Hanning  window  follows.  The 
subroutines  TRIMXD  and  EIGVLD  are  presented  in  reference  11,  and  sub- 
routines DPMCOS  and  DPMFFT  are  given  in  reference  12. 

In  order  to  execute  the  approximation  (41),  the  line  under  state- 
ment number  2 is  changed  to  CALL  PROBA(BT,  P,  Y) . This  subroutine  for 
the  Cubic  window  is  also  presented  below. 

IN  I LOLA  P 

UlMU'lSiO'.  ,<  ( b 1 ) .f  (rill  *2(200)  *YI,OPM<25) 

JA I A YI,0PiV-3. 7 1902, -3. 29053* -3, 09923 *-2.8761 6* -2, 57583,-2. 32635* 
4.-2. 1)5375 *-l,b4ifl5, -1.26155 *-.84162* -.52440*-. 2533b, 0,  *,25335*. 5244 
** .841 U2.1 ,2-i  156 * 1.64485 *2,05375* 2. 32635 *2.57563,2.87816* 3. 09023# 
>3.29053.3.71902/ 

0 = 1.4405825  i S HANNING 

CALL  •■lOOElbutt.f') 

call  bUHJtij  (2*0.  * YN0RM(1 ) * 3 , * YN0KM  ( 25 ) ) 

C AlL  OojCTut2*H50.i335.*2650.*2735.) 

CO  it  1£)T=3,0 
riT=2 , *» INT 

CALL  bLTSNo(,-:*39*2.  ) 

CAll  Ll;iEbu(2*0*0.  .YNORI'U) ) 

C An.  LINE  bu (2* l * 9. , YIIOR'V  (25)  ) 
call  LlNt.bG  ( 2 * 1 * 3 , * YMORF  ( 25 ) ) 

CALL  LlMt  auU*  1 * 3,  . VN0KP  (D) 

CALL  HUE.  ju  ( 2 * 1 * Q«  * YNOid'  ( 1 ) ) 

CALL  bl-.TSMG  (2*30*1.) 
uO  21  U=1*H 

call  LlnLbo't  2 * J . vJ*.,05*  Yli0K8(  l ) ) 

CAll  L INE  bo  ( / * 1 * J*  • 25  * YNOKM  2b ) ) 

DO  22  0=2  * 29 

CALL  UNtbo  ( 2 * 0 * 0 , * YN0kv  ( J ) ) 

’-■<  CALL  L INEbb  (2*1*3.*  YMOHV  ( J ) ) 

CAlL  SLTSI’Hi  ( i * 30  *2 . ) 

00  23  1=1,51 
Lb  a< 1)=.06*( I-t) 

ho  1 IPslM 
P=  (liT/C ) * IP 
SL=(I)T/C-1. ) /(H- 1) 

PRINT  c,  i)T  , iJ * SL 

2 POi, , -a  r (////•  B1  •: '613.8*  ' P ='14, » S/L  =*E13,8) 

CAll  PKot  Op  ( ,,T  ,P  * Y ) 

PRb.r  3,  t 

3 (•  UhmAT  ( /bb.2 :j  .8 ) 

DO  a 1=1 • bl 
Q=M1N<  T(I ) , ,999999) 

Q=MAX (9. .000001  ! 

4 YU)=TlNOHM(w**l) 

CALL  LlNtibG(2*fcl*X*Y) 

1 CONTINUE 

CALL  PAGED ( 2*0, 1*1) 

C-4 
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11 


2 

3 

1 


4 


5 


6 

7 


9 

a 


CONTINUE 
CAUL  EXITtiU1) 
END 


SUBROUTINE  pKOiiDP  (BT #P » ANS ) 

PARAMETER  MslOO  S MAXIMUM  NUMBER  OF  PIECES 
PARAMETER  N=2048*N4l=N/4+l 

DOUBLE  PRECISION  R(M»M}«Q(M)*B(M)*E(^)#w(M)»F(M)#GH(N)fGl(l<J»L0(N4 
»1>,C.ERROR»OO.V.P1»SL.TPE,XI2»PR.AT»OELXI.S»U,WXOXI 


INTEGER  P*P1 
DIMENSION  ANS < 1 » 

Cxi ,4405625600  9 HANNING 

IF(P.LE.M)  SO  TO  1 


JXM 

PRINT  2#  P'J 
FOHMAT(/»  P ='14»* 
00  3 J=lr5l 


IS  GREATER  THAN  V s' 13/) 


ANS(J)=-1. 

RETURN 

ERRQRsl.D-U 
DEC V=. 0600 

PI=3, 1415926535697932400 
P1=P-1 

SLs(BT/C-l,UO)/Pl 
DO  4 KSO » Pi 
D(K+l)sU<K+SL) 

00  5 Jsl*P 
00  5 KsifP 
L=ABS<J-K)+I 
R(J,K)=0<L) 

CALL  TRIMXU (P»M»R*U*B) 
CALL  EIGVcD<P»L>0»U,WrF> 

TPE=2.00/lPt*EHROH) 

XI2SI.D100 


PRsO.OO 

00  6 J=P»1,-1 

PR=PR*COG(L(U)  ) 

AT=1,UO/(P-.J  + 1.DO) 

S=P*EXP  ( AT*  ( LOG  ( TPE*4T I -PR ) ) 
XI2=MIN(X12#S> 

NFsN/2 


ULLXI=2»L»I)1*PI/«NF*ULLV) 

S=XI2/0ELXi 

JCs(S4l,00)/(jF 

NlsNF-1 

DO  6 K=0#"1 

5=0,00 

DO  9 JSOfJC 

SXS+WFIDXI ( (Kto*NF)*n£LXI ) 

GK(K*l)sS 

G1  (K  + UsO.UO 

CALL  OPMCOS(CO»NF) 

0=1,4427*1.05  (Nf-  ) + .5 

CALL  UPMFF  T (GRiGI  , CO  * U , -1 ) 

5=2 , DO*DELX I /P I 

L-'V  ID  W — i < . i 
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10  ANSi  t J)  = 1,DO-5*OH(J) 

IF<uF.EG.h)  RETURN 
UU  xi  J=l»46*5 

11  PRINT  12.  ANS(o)fANS(d+l)»AN5(J+2) »*MS(J+3),AnS(J+4) 

PRINT  12.  ANS<51) 

12  FORMAT (/5E20.8) 

NFsn 

60  TO  7 

FUNCTION  U(T>  f*  HANNING 
DOUBLE  PRECISION  T.S'i 

iFa.se. i.uo)  so  to  i 
SIx2,D0*PI*r 

US2.D0/3.U0*  (l.UO-T)*{l.DO+.r'''C*C05(SI))  + (.5tiO/Pl)*SIN(SI) 
RETURN 

1 U=0,D0 
RETURN 

FUNCTION  WFIUXKXJ 

DOUBLE  PRECISION  X.XTOP.AL.PE.BI » TEVP.Sx 

IF ( X «GT  ,0  • Uo ) 60  TO  1 

XTOPsl.DlOO 

WFIUXI=;SDO 

return 

X IF (X ,GT , XTuP)  00  TO  3 

ALsX.OO 
BE=-E(P)vX/p 

00  2 JXslM'X 
BIxE ( JI ) *X/P 
TEMP=Al+BE*,jI 
aEsBE-AL'-BI 

2 ALsTEMP 
SQsAL*AL+UE+UE 

IF  (SQ*(X*ERi<0R)**2«6T .10.00)  XTOPx^Ili  (XTOP.a  ) 
WFIUXI=-£IL/(S0AX) 

RETURN 

3 *Fiuxiso«uo 
RETURN 

ENO 

S.UBROUT INE  PROBA(UT.PiANS) 
double  PRECISION  C.O.Q.fclV.FllL 
INTEGER  F.Pl 
DIMENSION  ANS(l) 

CSX. 82009566  » CUBIC 

PlaP-1 

SL=(BT/C-X,)/P1 

B=X. 

00  1 Ksl.Pl 

X Bs0+2. * (X. -FLOAT (K)/P)*U(K*5L)**2 
CAPK=£,*P/E) 

PRINT  101.  CAPK 

101  FORMAT l / * CAPK  IS  'E15.8) 

BsP/B 

IBsB 

FBsR-lB 

CALL  GAMMA  ( 1 «+FD .G»$£.}2 ) 

SD=LOG(DBLE(G) ) 

UO  5 Ksi.lti 
DsFB+K 
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GOsGD+LOG(Q  ) 

00  3 K=1>5X 
V=,06*(K-1) 

IF( V,GT ,0» ) GO  TO  6 

ANS(K)=0* 

GO  TO  3 
BV=B*V 

ANS (K)sEXP (B*LOG (BV ) -B V+F  L1L ( DBLE(B+1 , ) , BV )-G0) 

CONTINUE 

RETURN 

PRINT  4 1 U 

FORMAT (/*  PROBLEM  AT  B = *E15,8> 

RETURN 

FUNCTION  FliL(A»XU) 

DOUBLE  PRECISION  SD*TO»AO*XD'*A 

SO=1,DO 

TDsl.DO 

AD= A- l, DO 

DO  1 K=1 » 1000 

TD=TD*XO/(AO+K) 

SDs&D+TD 

IF(ABS(TD) ,LE.1.D-8*ABS(S0) ) 00  TO  2 
PRINT  3, 

FORMATL/'  1000  TERMS*/) 

FllLsLOG(SD) 

RETURN 

FUNCTION  U{T)  « CUBIC 
IF  IT  »GE  , 1 . ) GO  TO  1 
U=1024./15l.*(l.-T)**7 
IF (T  ,GE,0,75)  RETURN 
U=U-fll92./15l.*<.75-T)**7 
IF (T  ,GE,0,5)  RETURN 
U=U+26672«/l51 » *( « 5-T ) ♦♦? 

IF (T .OE ,0,25)  RETURN 
U=U-57344./l51,*<.25-T)»"*7 

RETURN 

U=0. 

RETURN 

END 
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